Direct calculation of the solid-liquid Gibbs free energy difference 
in a single equilibrium simulation 
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Computing phase diagrams of model systems is an essential part of computational condensed 
matter physics. In this paper we discuss in detail the interface pinning (IP) method for calculation 
of the Gibbs free energy difference between a solid and a liquid. This is done in a single equilibrium 
simulation by applying a harmonic field that biases the system towards two-phase configurations. 
The Gibbs free energy difference between the phases is determined from the average force that the 
applied field exerts on the system. As a test system we study the Lennard- Jones model. It is shown 
that the coexistence line can be computed efficiently to a high precision when the IP method is 
combined with the Newton-Raphson method for finding roots. Statistical and systematic errors are 
investigated. Advantages and drawbacks of the IP method are discussed. The high pressure part of 
the temperature-density coexistence region is outlined by isomorphs. 



I. INTRODUCTION 

An important aspect of computational condensed mat- 
ter physics is to compute phase diagrams of model sys- 
tems. The naive approach is to preform a long-time sim- 
ulation at a selected state point and hope that the sys- 
tem by itself finds its preferred phase, i.e. the phase 
with the lowest Gibbs free energy. In most cases this 
strategy is not viable since first-order transitions are as- 
sociated with hysteresis. Thus the system is likely to be 
stuck in a metastable phase. The origin of this hystere- 
sis effect is the formation of an interface between two 
phases. The surface tension of the interface gives rise to 
a free energy barrier that the system has to overcome to 
transform from one phase to the other [Tj. A conceptual 
appealing and widely used strategy to overcome the hys- 
teresis problem is to preform simulations starting from 
an initial configuration with two phases in a periodic box 
[2HI1] • When a steady state situation is reached the sta- 
ble phase will grow at the expense of the other phase. 
The disadvantages of this approach are that: i) it relates 
to a non-equilibrium computation that cannot be done 
ad infinitum; ii) a sufficiently large system is needed to 
reach the steady state growth; iii) thermal fluctuation 
will result in some probability that the system will move 
towards the disfavored phase. In a recent paper 115^ these 
problems were resolved by introducing the so-called "in- 
terface pinning" (IP) method. In short, the idea of this 
method is to compute the average force needed to keep 
the system in the two-phases state. This is done by con- 
necting the system to a harmonic field which couples to 
an order-parameter that distinguished between the two 
phases of interest. The Gibbs free energy difference be- 
tween the phases is determined by the average force that 
the applied field exerts on the system. Thus, a standard 
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ad infinitum equilibrium simulation gives the information 
needed to computed the Gibbs free energy difference be- 
tween the two phases of interest. 

The purpose of this paper is to give a detailed descrip- 
tion of the IP method and show that it is a viable way 
of computing the solid-liquid Gibbs free energy and con- 
struct phase diagrams. As a test case, we investigate the 
Lennard- Jones (LJ) model [TB]. The remainder of the 
paper is organized as follows. In section II we describe 
the IP method in general terms. In section III we de- 
fine an order-parameter that distinguishes between solid 
and liquid by measuring long-range order. In section IV 
the IP method is applied to the LJ model, and statisti- 
cal and systematic errors are investigated. In section V 
we compare the IP method to other ways of computing 
Gibbs free energies and phase diagrams. The paper is 
completed with a summary. 



II. THE INTERFACE PINNING METHOD 

To introduce the IP method, imagine a two-phase sys- 
tem |17j in a periodic orthorhombic box elongated in the 
z-direction, i.e. with box lengths X <Y < Z (Fig. [I]). 
Consider configurations of the A^p^T-ensemble defined as 
the constant temperature and pressure ensemble where 
the box lengths X and Y are fixed at values so that the 
crystal is unstrained, while the box length Z is allowed to 
fluctuate in order to maintain a constant pressure. Two 
interfaces will form orthogonal to the long axis. This 
orientation will minimize the interface area and thereby 
the interface Gibbs free energy G^. The system is only 
barostated in the z-direction since the surface tension 
will add an a priori unknown pressure to the Px and Py 
components of the pressure tensor [18j. (It is worth not- 
ing that an orthorhombic box is not a requirement. The 
angle between the box vectors X and Y may differ from 
90°, but should then be kept constant at an angle not 
does not strain the crystal). 
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FIG. 1. Two-phase configuration of the LJ model in a 
periodic orthorhombic box at a state point where the liq- 
uid is the thermodynamically stable phase while the crystal 
is metastable. This is an equilibrium configuration since a 
harmonic field biasing towards two-phase configurations have 
been applied. The average force exerted by the applied field 
on the system relates to the Gibbs free energy difference be- 
tween the phases. 



Assume that the system is sufficiently large so that 
the central regions of the pure phase slabs exhibit bulk 
properties up to some arbitrary threshold values. Parti- 
cles can then be labeled either s = [solid], I = [liquid] or 
i = [interface] and the total number of particles can be 
written as 

N^Ns+Ni + N,. (1) 

Let /Xs and fii be the chemical potential of the solid and 
the liquid, respectively. The total Gibbs free energy of 
the two-phase system is then 

G^Nsfis + Ni^ii + G^. (2) 

as sketched on Fig. [2] When the relative position of 
the interfaces change within the two-phase regime, Gi 
and Ni can be replaced by an constant if "wiggles" can 
be neglected (wiggles [19l [20] are discussed later in the 
paper). Thus combining the last two equations gives 

G = NsAfi + const. (3) 

where A/i = fig — fii. Throughout the paper we let "A" 
denote "[solid] — [liquid]" and let "const." refer to an 
constant. 



A. Harmonic field biasing towards two-phase 
configurations 

To maintain the system in configurations having two 
phases, i.e. "pinning the interfaces" , we apply a harmonic 
field that couples to an order-parameter which relates to 
the amount of crystal particles in the simulation box. To 
this aim we introduce a global order parameter Q(R.) 
where R = {ri, r2, . . . , r^} is the configuration of par- 
ticles (for simplicity we assume that particles are point- 
like). Let Qs and Qi be the average values of (5(R) when 
the system is completely solid or liquid, respectively (at 




liquid two phases crystal 

Order-parameter Q 



FIG. 2. Sketch of Gibbs free energy G{Q) (solid; black) along 
an order-parameter Q that measure the amount of crystalline 
particles in a system similar to the on shown on Fig. [T] The 
liquid have the lowest G and it is thus the thermodynam- 
ically stable phase while the crystal is a metastable phase. 
The double arrows in the center of the figure indicate the in- 
terface contribution Gi to G{Q). The double arrow on the 
right hand side of the figure indicates the total change of G 
when moving from one phase to the other. The dashed curve 
indicates the Gibbs free energy G'{Q) of a system where a har- 
monic external field has been applied to stabilize two-phase 
configurations. The inset shows a computed G{Q) for the LJ 
model in the two-phase region (A'' — 5120; T = 0.8; p = 1.5; 
computed using umbrella sampling 21] )■ 



a given pressure p and temperature T). We define Q so 
that it has a linear dependency on the amount of solid 
particles in a two-phase state when additional degrees of 
freedom are integrated out (such as phonon vibrations in 
the slabs of the pure phases): 

Q - ^Qs + ^Qi + = + const. (4) 

where is a constant contribution from interface par- 
ticles. Let ?7(R) be the energy of the unperturbed sys- 
tem, and 

t/'(R) = [/(R) + |[g(R)-a]2, (5) 

be the energy when a "spring-like" harmonic field is ap- 
plied. We refer to the field parameters a and k as the 
anchor point and the spring constant, respectively. The 
Gibbs free energy along the Q coordinate when the bias- 
ing field applied is (dashed line on Fig. [2]) 

G'(Q) = G(Q) + ^[Q-a]2. (6) 

To give an expression for the probability distribution of 
Q we use that P'{Q) oc exp{—G'{Q)/kBT). By insertion 
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FIG. 3. P'iQ) distribution of a two-phase system where the 
relative interface position have been pinned with a harmonic 
field (LJ model; T = 0.8; p = 1.5; = 2.5; tsim = 4000). Four 
values of spring constants have been used, k = {2, 4, 10, 20} 
respectively. Fluctuations of the order-parameter Q follow 
Gaussian statistics (dashed lines; Eq. ([7|). The liquid is the 
thermodynamically stable phase at this state point and the 
average value of Q is pulled by the system to values below the 
anchor point of a = 27. 



of Eq. ([3]) and elimination of Ng with Eq. Q we get 
that. 



Q - a 



kAQ 



(7) 

in the two-phase regime. This distribution is shown for 
the LJ model in Fig. |3] 



C. Computing coexistence state points with the 
Newton-Raphson method for finding roots 

Coexistence state points are defined as Afi{p,T) — 
and may be computed efficiently using the Newton- 
Raphson algorithm for finding roots. The required 
derivatives of A/i along isobars and isotherms are given 
by the standard thermodynamic expressions 



dm 



and 



where 



dp 



dT 



Av 



-As 



As = 



Au + pAv — A^ 



(10) 



(11) 



(12) 



In these relations. Ad, As, and Au are changes in specific 
volume, entropy, and energy, respectively. Thus coexis- 
tence points can be computed by iteration of 



or 



A^' 



Ai;(*) 

A^(') 
As(^) ■ 



(13) 



(14) 



Iterations are continued until the computed A/i is zero 
within the statistical error. 



D. Algorithm for computing coexistence state 
points 



B. Computing Afj, from the average force exerted 
by the applied field on the system 

The chemical potential difference A/x can be computed 
from the average force 

= -«[(Q)' - a] (8) 

that the field exerts on the system (along the Q coordi- 
nate). When equilibrium is established the relative po- 
sition does not change up to thermal fluctuations and 
^ficid ^ _^systom ^^^^.^ ^system = - |g . By applying 

the chain rule Afi =-§§- = ^-^^ then 

Am = -^[(Q)'-«] (9) 

where = is obtained from Eq. Alternatively, 
a statistical mechanical deviation of tms is possible by 
isolating A/i from the average of the P'{Q) distribution 
(Eq. 0). 



To conclude this section, we give an algorithm for com- 
puting coexistence state points in the phase diagram: 
First, select a pressure p and temperature T for the initial 
set of simulations. Then, 

i: Construct a crystal configuration in an elongated or- 
thorhombic box; 

ii: Determine the lattice constants of the unstrained crys- 
tal by preforming NpT simulation where the box 
lengths X, Y and Z are allowed to fluctuate inde- 
pendently to maintain constant pressure; 

iii: Compute Qg and Vs in an Np^T simulation of the 
unstrained crystal; 

iv: Construct a liquid configuration in an elongated or- 
thorhombic box having the same box lengths X and 
Y as the unstrained crystal; 

v: Compute Qi and vi in a Np^T simulation of the liquid; 
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vi: Construct a two-phase configuration having the same 
box lengths X and Y as the unstrained crystal; 

vii: Compute (Q) in an Np^T simulation of the two- 
phase system with an interface pinning |((5(R) — 
a)^ field applied; 

viii: Calculate A/i using Eq. (|9]); 

ix: If is non-zero within the statistical error, repeat 
steps i-ix at the pressure given by Eq. (131 or the 
temperature given by Eq. ( 14 ) . 



We note that an algorithm only involving two-phase sim- 
ulations can be designed, since a two-phase simulation 
contains information about the bulk properties of the liq- 
uid and the crystal. In this paper we choose to use the 
above algorithm for practical reasons. 



III. TRANSLATIONAL ORDER PARAMETER 

To utilize the method we need to define an order pa- 
rameter Q(R) that distinguishes between the phases of 
interest. Unlike liquids, crystals have long-ranged trans- 
lational order and the collective density field may be used 
to define (3(R): 



Q(R) = |Pk| 



where 



N 



|exp(-ik- rj), 



(15) 



(16) 



where Fj is the force without external field, and 

, I ,^ jR[pk] sin(k ■ rj) + S[pk] cos(k ■ rj) 

VjiPkl - -k — (18j 

\pk\VN 



where 5i[pk] = X^jLi cos(k • Vj)/^/N and 



5[Pk] = 



— X]j=i sin(k • Yj)/yN are the real and imaginary parts 
of pk respectively. Forces can be computed with an effi- 
cient 0{N) scaling algorithm although the force on par- 
ticle j depends on the position of all of the particles (this 
typically result in a 0(7V^) scaling algorithm). This is 
done in two 0{N) steps: i) compute pk using Eq. (16) 



and ii) compute particle forces using Eqs. (171 and (18 1 



Monte Carlo sampling involves evaluation of the energy 
change SU' when a particle is moved or the box length 
Z is changed: 



SU' = SU + -SIp^]-" ~ KaS\pk\ 



(19) 



where S\pi^\ 



|Pk 



try 12 



I ^current 1 2 



|Pk 



try I 



I ^current I ^ These chaugcs may be computed by evaluating 
Sp^ = p^^y - pcurront ^j^g current value of pk = p™'™"* 
is stored. Moving particle j yields 



Spy, = [exp(-ik • rj^) - exp(-ik • r™"™*)]/V7V. (20) 

Thus computing Spk only involves information about par- 
ticle j allowing for efficient computations. When the box 
length Z is change Spk = and SU' = SU since k is 
perpendicular to the z-direction. 



and k = {2TTnx/X,2'Kny/Y,0). The integers and Uy 
are chosen so that the wave vector k corresponds to a 
Bragg peak. This will maximize the contrast between 
the liquid and the crystal. The k vector is in the xy- 
plane (n^ = 0) since Z fluctuates in the A'^p^T-ensemble. 
The factor N~2 ensures scale invariance for the average 
liquid value, Qi oc 1, while the intensity of a (single) 
crystal will scale as Qg oc N^. For two-phase configu- 
rations Q will scale linearly with the amount of crystal 
particles (fulfilling Eq. ^) since the pk-argument of the 
crystal slab and liquid slabs is independent of each other. 
We note that this order-parameter may be problematic 
in the supercooled regime, since a crystal can lower |pk| 
by introducing a long-wave length displacement of par- 
ticles. This can be avoided by using a different order- 
parameter, e.g. the Steinhard Qq order-parameter |22j . 
This was done in Ref. We choose to use |pk| as 

order-parameter since it is conceptually appealing, sim- 
ple and general applicable. 

Equilibrium trajectories can be constructed using stan- 
dard Monte Carlo sampling or standard Molecular Dy- 
namics simulations [21] . For the latter, forces exerted on 
particles by the external field have to be evaluated: The 
force acting on particle j is 



F' 



K(|Pk| - a)Vj|pk| 



(17) 



IV. SOLID-LIQUID COMPUTATIONS OF THE 
LENNARD-JONES MODEL 

As a test case we apply the IP method to compute 
solid-liquid Gibbs free energy differences of the LJ model 
LJ interactions are truncated and shifted: U = 

for r.ij < Vc and zero otherwise. LJ units are used 
throughout the paper: e = a = m = kB = l. Two 
truncation distances are considered: Tc = 2.5 and Tc — 6. 
The first choice of 2.5 is the standard truncation of the LJ 
model. The latter choice of 6 is a better approximation 
of the full LJ model (r^ — > oo). Molecular dynamics sim- 
ulations are perform using the LAMMPS software pack- 
age [53]. The pk-field was implemented into the package. 
The Parrinello-Rahman barostat is used [24] with a time 
constant of 8 Lcnnard- Jones time units together with a 
Nose-Hoover ^25^ 23^ thermostat with a time constant of 
tnh — 4. Trajectories are evaluated using a time step of 
0.004. 

As an example, we compute Ap at p — 1.5 and T — 0.8 
(tc — 2.5) as described in the following. First, a crys- 
tal structure of 8x8x20 face centered cubic unit cells 
{N ~ 5120) is constructed and simulated for tgim = 4000. 
All box lengths are allowed to fluctuate in order to de- 
termine the lattice constants of the unstrained crystal 
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FIG. 4. Panels (a) and (c) show the Gibbs free energy dif- 
ference between solid and liquid computed with IP method 
(+) and thermodynamic integration (lines). Panels (b) and 
(d) show the specific entropy and the specific volume, respec- 
tively. The lines on the lower panels are cubic polynomial fits. 
The lines in the upper panels are computed by integration of 
the fits. The integration constant is chosen to give the best 
overall agreement {r^ = 2.5). 



FIG. 5. Coexistence region of the Lennard-Jones model in 
the pT-plane on a log-10 scale. Filled symbols are computed 
with the IP method (Tables [T] and T he points labeled as 
X 's are reproduced from Reference [27] . The agreement is 
good. At high temperatures and densities, the coexistence 
region is outlined by isomorphs (see text for details). The 
shape of the isomorphs are determined at the state points 
indicated by open circles (no fitting procedure was applied). 



giving box lengths of X = Y = 12.92. The unstrained 
crystal is then simulated for tgim — 4000 in the Np^T 
ensemble, and Qs = 55.04 {ux = 16, Uy = 0) and the 
average partial volume Vg — 1.052 is recorded. Next, a 
liquid configuration is constructed by melting the crys- 
tal in a constant volume simulation by simulating at a 
high temperature (T = 5). The iVp^T-ensemble (us- 
ing X = Y = 12.92) of the liquid is then simulated for 
^sim = 4000. Qi = 0.93 and an average specific volume 
of vi — X.YIl is recorded. Then, a two-phase configu- 
ration is constructed by performing a high temperature 
constant volume simulation where particles at z < Z/2 
are kept at their crystal positions using harmonic springs 
anchored at crystal sites. The box volume is set in be- 
tween that of the crystal and the liquid by scaling the box 
length Z . The iVp^T-ensemble with a harmonic bias-field 
(a = 27; k = 10) is then simulated for fgim = 40000 to 
compute (Q)' — 25.246. Eq. ([9]) yields a chemical po- 
tential difference of A/i — 0.080. A configuration from 
this last simulation is shown in Fig. [T] The two up- 
per panels in Fig. |4]show A/i along the p = 1.5 isobar 
and the T — 0.8 isotherm, respectively [tc = 2.5). The 
solid lines are computed with thermodynamic integration 
of As and Av, respectively (shown in the lower panels). 
The agreement is excellent. 

Next, we use the Newton-Raphson method along 
isotherms to compute coexistence state points. As an 
example we computed the T = 0.8 co exis tence pressure 
from the state point at p^^^ = 1.5. Eq. 13 provides pres- 
sures of = {2.141,2.189,2.185(2)}. In the last iter- 



ation the estimated chemical potential difference is zero 
within the statistical error, A/j, — 0.0001(2) (numbers 
in parentheses indicate the statistical errors on the last 
digit) . Table |l] and |ll] list computed coexistence points 
using Tc = 6 and = 2.5, respectively. As a consistency 
check, we note that the computed melting line obey the 
Clausius-Clapeyron relation, — ^ (two last rows 
in Tables |l] and |ll| . The left-hand side of the relation is 
computed by central differences of the computed melting 
line. 



A. Isomorph prediction of the pT coexistence 
region 

an aside, we test a recent theoretical prediction 
related to the melting region in the pT-plane (Fig. 
I5|: a large class of systems have curves in the phase di- 
agram, referred to as "isomorphs" [28], where structure, 
dynamics and some thermodynamic properties are nearly 
constant. These are defined as 




(21) 



where T* is the temperature at a reference state point 
and h(p) [3D] is a function of p (not to be confused with 
the specific enthalpy h). This class of "simple" [3T] sys- 
tems are characterized by the property that virial W (the 
potential part of pressure) and potential energy U fluc- 
tuations are strongly correlated in the NVT ensemble 
ifSW = W- {W) and SU = U - (U) then the 
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TABLE I. Solid-liquid coexistence line of the truncated LJ model {rc — 2.5) 



Ttn 


0.6 


0.7 


0.8 


0.9 


1.0 


1.2 


1.4 


1.6 


1.8 


2.0 


2.2 


2.4 


2.6 


Pm 


-0.212 


0.928 


2.185 


3.514 


4.939 


7.921 


11.181 


14.632 


18.180 


22.007 


26.029 


30.050 


34.314 




-1.046 


0.049 


1.264 


2.555 


3.943 


6.859 


10.056 


13.448 


16.943 


20.717 


24.688 


28.661 


32.878 


Vs 


1.0614 


1.0452 


1.0277 


1.0110 


0.9951 


0.9672 


0.9421 


0.9202 


0.9014 


0.8835 


0.8671 


0.8530 


0.8394 


1)1 


1.2194 


1.1714 


1.1360 


1.1080 


1.0830 


1.0446 


1.0117 


9838 


0.9612 


9399 


0.9211 


0.9043 


8888 


Ua 


-5.358 


-5.156 


-4.953 


-4.742 


-4.513 


-4.020 


-3.483 


-2.907 


-2.301 


-1.663 


-0.997 


-0.315 


0.394 


Ul 


-4.294 


-4.218 


-4.075 


-3.888 


-3.683 


-3.183 


-2.627 


-2.041 


-1.400 


-0.727 


-0.009 


0.696 


1.446 


As 


-1.718 


-1.507 


-1.392 


-1.327 


-1.263 


-1.207 


-1.168 


-1.123 


-1.107 


-1.091 


-1.087 


-1.063 


-1.055 


dTm 




12.0" 


12.9 


13.8 


14.7 


15.6 


16.8 


17.5 


18.4 


19.6 


20.1 


20.7 




As 
Av 


10.9 


11.9 


12.9 


13.7 


14.4 


15.6 


16.8 


17.6 


18.5 


19.3 


20.1 


20.7 


21.4 



computed by central difference of values in tlie first two rows. 



TABLE 11. Solid-liquid coexistence line of tfie truncated LJ model (r^ = 6) 





0.6 


0.7 


0.8 


0.9 


1.0 


1.2 


1.4 


1.6 


1.8 


2.0 


2.2 


2.4 


2.6 


Pm 


-0.970 


0.132 


1.337 


2.629 


4.012 


6.930 


10.145 


13.549 


17.104 


20.857 


24.850 


28.916 


33.041 


Pm+P^""'^ 


-1.030 


0.068 


1.270 


2.560 


3.940 


6.853 


10.063 


13.463 


17.014 


20.763 


24.753 


28.815 


32.937 


Vs 


1.0553 


1.0400 


1.0242 


1.0086 


0.9933 


0.9661 


0.9412 


0.9194 


0.9002 


0.8827 


0.8663 


0.8518 


0.8388 


Vl 


1.2358 


1.1804 


1.1425 


1.1120 


1.0864 


1.0470 


1.0127 


0.9852 


0.9616 


0.9403 


0.9208 


0.9038 


0.8885 


Us 


-6.314 


-6.125 


-5.929 


-5.722 


-5.506 


-5.037 


-4.526 


-3.972 


-3.385 


-2.768 


-2.120 


-1.451 


-0.764 


Ul 


-5.024 


-5.008 


-4.902 


-4.755 


-4.570 


-4.106 


-3.603 


-3.021 


-2.409 


-1.762 


-1.085 


-0.379 


0.325 


As 


-1.858 


-1.622 


-1.480 


-1.377 


-1.311 


-1.245 


-1.177 


-1.153 


-1.125 


-1.103 


-1.085 


-1.073 


-1.050 


dpm, 
dT^ 




11.5" 


12.5 


13.4 


14.3 


15.3 


16.5 


17.4 


18.3 


19.4 


20.1 


20.5 




As 
Av 


10.3 


11.6 


12.5 


13.3 


14.1 


15.4 


16.5 


17.5 


18.3 


19.1 


20.0 


20.6 


21.1 



dT^ computed by central difference of values in tfie first two rows. 



correlation coefBcient R ^ {SW5U} /yW^VfWW) is 
close to unity. 

For systems with inverse power-law pair interactions, 
Uij oc r~", isomorph invariance is trivial |35j with h(p) = 
(p/p*)"/'^ where is the density at the reference state 
point. For systems with LJ pair interactions (a sum of 
two inverse power-laws) the scaling is approximate and 
reflects an effective inverse power-law of repulsive inter- 
actions |36) . The apparent exponent is state point and 
phase dependent [301 EZ] : 



h(p) = 



- 1 



- 2 



(22) 



where 7* is a constant that may be determined from 
virial-energy fluctuations in the NVT ensemble at the ref- 
erence state point by using the identity 7 — q'°|(p) = 

^{{suy) • Here it is used that excess entropy Sex is iso- 
morph invariant (7* = 7 evaluated at the reference state 
point). 

The dashed line on Fig. [s] is a liquid isomorph (7* = 
4.816; = 2; p* = 1/0.9403 = 1.064; R = 0.991) and 
the solid line is a crystal isomorph (7* — 5.517; = 2; 

= 1/0.8827 = 1.133; R = 0.998). At high temper- 
atures and densities the coexistence region is outlined 



by these isomorphs. Deviations from the prediction are, 
however, significant at low pT. This is properly due to 
long-range attractive interactions not being accounted for 
by an effective inverse power-law. Consistent with this 
interpretation, the melting region of the = 2.5 trun- 
cation of pair interactions deviates from the rc = 6 in 
this part of phase space. A scaling form of Ap^ + Bp^ 
(like Eq. (22)) has previously been validated by Khra- 
pak and coworkers [351 ■ They used Rosenfeld's rule of 
additivity of melting curves |40| to motivate the scaling 
law. 



B. Correcting for missing long-range attractions 

To estimate the melting line of the full LJ model 
— ^ 00), we applying an approximate pressure cor- 
rection pi-^il that rectifies missing long-range attractions 
of the truncated model. To this aim we first consider the 
pressure correction in a simulation of solid or liquid in 
bulk: it is convenient to assume that the radial distri- 
bution function is constant at distance larger than the 
truncation. Then the correction is analytic and only de- 
pends on p and Tc [2T]. Since the densities of the solid 
and the liquid are different, so are the corrections for the 
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FIG. 6. Coexistence line of the Lennard- Jones model in the 
pT-plane. Filled symbols are coexistence points computed 
with the IP method (Tables |l] and |lT| corrected for truncation 
of long-range contributions to the pressure, Eq. ( |23[ ). The 
solid line is a cubic fit (rc = 6). The symbols +'a and x's are 
coexistence points reproduced from References [JT] and [27j . 
respectively. The inset shows deviations from the fit. The 
asterisk is the gas-liquid critical point (Tqp ~ 1.31; p^p = 
0.15) f?51. 



FIG. 7. Q time autocorrelation function for four spring 
constants k (same parameters as results shown on Fig. [sf. 
Decorrelation occurs on two timescales: i) a fast time scale 
related to sound waves and ii) interface movements. Dashed 
lines are Aexp(— t/r) fits to the slow interface process. The 
inset shows the relative statistical error on the = 0.080 
estimate. This error is computed by dividing runs into inde- 
pendent blocks t43j of length tbiock = 100 (4im ~ 2000). 



two phases. For the pressure correction of a two-phase 
simulation we use the average of the bulk pressure cor- 
rections: 



„tail 



Svr 



— r 
3 " 



(23) 



Statistical error 



How does the statistical error of the A/i estimate de- 
pend on the choice of spring constant k7 To answer 
this question, we compute the Q{t) autocorrelation func- 
tion (using the Wiener-Khinchin theorem ^44^ with Fast 
Fourier Transforms) 



Third rows in Tables |I] and |IT] list the corrected melting 

pressures {pm +p^^^^). Deviations between the corrected 
melting pressures when truncating at r^, = 2.5 or = 6 
are comparable to statistical error (Tables |l] and |ll|. 

Computed melting points are shown on Fig. [6] as filled 
symbols. In the same figure, +'s and x's are coexistence 
points computed with other methods [271 HI] • The agree- 
ment is excellent. Differences are highlighted in the insert 
by showing deviations from a cubic fit to the computed 
melting points. Deviations from the results of Reference 
[27j are within statistical error, while the melting pres- 
sure of reference [H] is systematically to low by about 
Ap ~ 0.05 (except for one data point). Systematic er- 
rors in computed melting lines are common [IT] and are 
typically related to approximate tail corrections (like Eq. 
(23)), finite size effects or method specific systematic er- 
rors. In the following subsections we will discuss system- 
atic and statistical errors related to the IP method. 



C{t) 



{5Q{0)5Q{t)) 

my) 



(24) 



where 6Q{t) = Q{t) - (Q). Fig. [t] shows C{t) for four 
choices of k (p = 1.5; T = 0.8; = 2.5). C{t) re- 
veals two relaxation processes that occurring on different 
time-scales. We assign them as follows: i) a fast pro- 
cess related to phonon vibrations and rearrangements of 
particles, and ii) a slower over-damped process related to 
particles moving between phases, Ns{t) = —Ni{t)+cox\sl. 
For the investigated k's the characteristic time r for the 
slow process is nearly constant, r scales as 1/k for smaller 
values of k this timescale (data not shown). The rela- 
tive statistical error on the computed is estimated 
by dividing runs into 20 independent blocks of length 
^biock = 100 > r |43j . For the shown k's, spanning two 
orders of magnitude, the statistical error is independent 
of K (inset on Fig. [7|. 
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FIG. 9. Two two-phase configurations of a square lattice gas 
with attractive interactions between neighbor particles. Solid 
particles are colored dark gray (red) and fluid particles are 
colored light gray (green). The interface Gibbs free energy d 
is different for the two configurations and moving a particles 
from one phase to the other, as indicated by the arrow, result 
in different changes of the energy. This result in wiggles on 
G{Ns). 



FIG. 8. Computed Ajj. (Eq. [9| as a function of the "interface 
distance" defined as ((Q)' - Qi)/[Qs - Qi) [T = 0.8; p = 1.5; 
Tc = 2.5; K — 10) of three system sizes with the same geom- 
etry ({4x4xl0;6x6xl5;8x8x20}; iV = {640,2160,5120}). 
The inset shows the average of the computed Afi's as a func- 
tion of system size. Error bars indicates the statistical error 
plus the variation related to interface position (tsim = 40000). 



D. Systematic errors at small system sizes 

To investigate finite size effects we computed A/z at 
p — 1.5 and T = 0.8 with Vc = 2.5 using three systems 
sizes with the same geometry: N — {640, 2160, 5120} re- 
spectively (Fig. [8|. For the smallest system sizes the 
computed A/i depends measurably on the relative inter- 
face positions (varied by changing a). This effect is not 
seen for the two larger system sizes. The inset shows 
the system size dependency of the compute A/i. Error 
bars indicate statistical error plus the variation related to 
the positions of the interfaces relative to each other. The 
dashed 1 /N line is a guide to the eye. For the largest sys- 
tem sizes, the error on the estimated Afi is on the order 
of a 10~^. This correspond to an error o n th e computed 
melting temperature of about lO^'^ (Eq. 14 1. 



E. 



Gibbs free energy dependency of the interface 
positions relative to each other 



We have assumed that the Gibbs free energy in the 
two-phase region is independent of interface positions rel- 
ative to each other. There are, however, two effects that 
may spoil this assumption: i) if the distance between the 
interfaces is sufficiently small, particles in one (or both) 
phases will not have bulk properties, and ii) "wiggles" on 
G{Ns) dH |20]. To exemphfy the latter effect think of a 
square lattice gas with attractions between neighboring 
particles. Fig. |9] shows two two-phase configurations of 
this model. The interface Gibbs free energy Gi is dif- 




FIG. 10. Crosses show G{Q) in units of fcsT of six ap- 
plied fields with different spring anchor points (data points 
are shifted horizontally for clarity). The circles show G{Q) 
computed with umbrella sampling ^1]. G{Q) has a linear 
dependency on Q (within the statistical noise). 



ferent for the two configurations resulting in wiggles of 
G{Ns). The wriggle period on G{Q) is AQ/N^ where 
Nz is the number of crystal planes in the ^-direction. To 
investigate the Gibbs free energy dependency of interface 
positions of the LJ model, we perform simulations over a 
range of a's with overlapping P'{Q) distributions. From 
this G{Q) is constructed with the umbrella method [2T] 
(histogram reweighing is done with the MBAR algorithm 
[45]). We find no wiggles (Fig. 



10) 



but conjecture that 
they are hidden in the statistical noise. We emphasize 
that wiggles may be accounted for and do not constitute 
a fundamental limitation of the IP method. 



F. Guidelines for choosing a and k 

How should the anchor point a and the spring constant 
K of the harmonic field be chosen to yield the optimal 
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computation? To answer this question, we note that the 
average distance between the two interfaces should be 
large as possible to ensure that slabs of the pure phases 
have bulk properties. This distance can be controlled by 
the anchor point a. To give a guideline for an optimal a, 
we consider the limit where k — >■ oo or A// = 0. From 
Eq. [9] we find that a — (Q)' . The optimal value of a is 

a^Qi + ^. (25) 

The "~" indicates that the interface contribution Qi has 
been ignored. Next, we consider the spring constant k. 
There are two arguments for choosing a stiff spring, i.e. 
large value of k: i) a small k would not keep the system 
in two-phase configurations; ii) the relaxation time of in- 
terface dynamics r scales as 1/k for small k's, thus giving 
bad statistics. There are, however, also an argument for 
choosing a small k: Interface fluctuations should span at 
least one crystal plane to account for wiggles on G{Q). 
Thus, K should optimally be chosen to that interface fluc- 
tuation span one crystal plane. Setting the standard de- 
viation of the P'{Q) distribution (Eq. equal to AQ/N^ 
we get: 

- - kBT^, (26) 

where Nz is the number of crystal planes in the z- 
direction. We emphasize that one of the conclusions of 
this paper is that the IP method is forgiving towards the 
choice of held parameters. 

V. ADVANTAGES AND DRAWBACKS OF THE 
INTERFACE PINNING METHOD 

Let us briefly review other methods for computing 
Gibbs free energies and phase diagrams [HI [21] be- 
fore discussing the advantages and drawbacks of the IP 
method. In the moving interface approaches, discussed in 
the introduction of the paper, a simulation is performed 
of a two-phase system pt414|. The thermodynamical fa- 
vored phase will grow in a constant NpT or ^VT simu- 
lation, allowing to locate coexistent points by changing 
intensive variables such as p, T or fx. An alternative is 
to use an indirect method where the Gibbs free energy 
of the pure phases is computed in separate simulations. 
The Gibbs free energy can be computed by Widom inser- 
tion [46l Elj or by thermodynamic integration to a state 
of know Gibbs free energy, e.g. an ideal gas [3S], a har- 
monic solid [5S] or an Einstein solid [49j. Umbrella sam- 
pling or metadynamics along a good reaction coordinate 
can be used to compute the Gibbs free-energy of trans- 
forming the system from one phase to the other 50 -.53J. 
The reader is encouraged to explore Refs. [TTl [5T] and 
references within for more about methods for computing 
phase diagrams. 

What are the advances and drawbacks of the IP 
method? First, the IP method inherits the conceptual 



simplicity, general applicability and ease to implement 
of the moving interface method. Since the solid-liquid 
interface is represented explicitly, simulations give infor- 
mation about interface properties. The surface tension 
may be computed by integration of pressure tensor el- 
ements [18] or from capillary fluctuation [M]. Crystal 
growth rates may be computed from Q{t) fluctuations 
(Fig. [7]) similar to the method suggested by Briels and 
Tepper [SS] ES] (in the BT method, two-phase configura- 
tions are stabilized by keeping the volume constant). We 
leave such investigations to future publications. A dis- 
advantage of having an explicit interface is that it is in 
direct contact with the solid and the liquid phases [57] . 
In effect, properties of the liquid and the solid slabs may 
not have bulk values. This is evident in the computed 
A/x's of the system size with 640 particles (Fig. [s]). Here 
A/x depend systematically on a. This may lead to larger 
finite size effects compared to methods where free ener- 
gies are computed in separate simulations having periodic 
boundaries. Another disadvantage is that a low interface 
mobility can result in slow dynamics (Fig. [7|). This will 
result in large statistical errors or, in the worst case, it 
may even be difficult to reach equilibrium. 



The IP method constitutes an alternative when other 
methods are difficult or impossible to use. Widom in- 
sertion 137] works best for dilute systems whereas it 
is not a viable option for dense liquids or solids. Ther- 
modynamic integration to a state of known Gibbs free 
energy El] is problematic when the path include addi- 
tional first order transitions. This happens when a phase 
is surrounded by other phases in phase diagram. A refer- 
ence state point can also be nontrivial to identify. Exam- 
ples are quasi-crystals, liquid-crystals, plastic-crystals or 
other phases having a mixture of order and disorder. The 
IP method is versatile, and may be used to study these 
phases. Moreover, it can be generalized to be used with 
two-phase simulations of only fluid phases (gas-liquid or 
liquid-liquid) or only solid phases. For the latter, one of 
the crystals will be strained when simulating the NpzT 
ensemble of a two-phase configuration (if the lattice con- 
stants of the crystals are different). This can be rectifled 
by adding the free energy of straining to the computed 
A/i. Low mobility of solid-solid interfaces could make it 
unfeasible to use the IP method. Integration along a re- 
action path of transformation is not trivial, since it often 
is difficult to identify a suitable coordinate capturing the 
entire phase transition. As an example, the number of 
crystalline particles Ns is not a good reaction coordinate, 
since the cluster of crystalline particles undergoes a ge- 
ometric transition from spherical when Ns is small to a 
slab when ~ N/2 (Fig. [T]). Thus, there is at least one 
geometrical free energy barrier orthogonal to the Ng co- 
ordinate (with a barrier height that scales with the area 
of the cluster: N^ ). With the IP method, the selection 
of the order parameter Q only has to distinguish between 
the two phases of interest. 
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VI. SUMMARY 



In summary, we have given a detailed description of 
the IP method and shown that it can be used for effi- 
cient calculations of the Gibbs free energy difference be- 
tween a solid and a liquid. The melting line can be com- 
puted efficiently to a high precision when the method is 
combined with the Newton-Raphson algorithm for find- 
ing roots. As an example, the solid-liquid coexistence 
line of the truncated LJ model line was computed. As 
an aside, it was shown that the high pressure part of the 
temperature-density coexistence region is outlined by iso- 
morphs. An approximate pressure correction for rectify- 
ing truncation of pair interactions was given. Statistical 
errors and systematic variations were investigated. 

An important advantage of the IP method is that the 



solid-liquid Gibbs free energy difference is computed di- 
rectly in a ad infinitum simulation at a single state point. 
This makes it versatile and a viable alternative when it is 
difficult or impossible to preform Widom insertion, ther- 
modynamic integration or integration along a reaction 
path of transformation. 
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